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Abstract 

Historical and contemporary evolutionary processes can both contribute to pat- 
terns of phenotypic variation among populations of a species. Recent studies are 
revealing how interactions between historical and contemporary processes bet- 
ter explain observed patterns of phenotypic divergence than either process alone. 
Here, we investigate the roles of evolutionary history and adaptation to current 
environmental conditions in structuring phenotypic variation among polyphenic 
populations of sunfish inhabiting 12 postglacial lakes in eastern North America. 
The pumpkinseed sunfish polyphenism includes sympatric ecomorphs specialized 
for littoral or pelagic lake habitats. First, we use population genetic methods to 
test the evolutionary independence of within-lake phenotypic divergences of eco- 
morphs and to describe patterns of genetic structure among lake populations that 
clustered into three geographical groupings. We then used multivariate analysis 
of covariance (MANCOVA) to partition body shape variation (quantified with 
geometric morphometries) among the effects of evolutionary history (reflecting 
phenotypic variation among genetic clusters), the shared phenotypic response of all 
populations to alternate habitats within lakes (reflecting adaptation to contempo- 
rary conditions), and unique phenotypic responses to habitats within lakes nested 
within genetic clusters. All effects had a significant influence on body form, but 
the effects of history and the interaction between history and contemporary habitat 
were larger than contemporary processes in structuring phenotypic variation. This 
highlights how divergence can be better understood against a known backdrop of 
evolutionary history. 



Introduction 

Trait differences among populations and species are often 
interpreted as being the result of contemporary evolutionary 
processes, such as adaptation to current ecological condi- 
tions, or the result of historical evolutionary process such 
as bottle necks, genetic drift, or adaptation to past ecologi- 
cal conditions. Comparative approaches are frequently used 
to test for adaptation to ecological conditions by demon- 
strating a significant correlation between current conditions 
and organismal traits in an attempt to reject random chance 
as an explanation for trait differences among populations 
(Robinson et al. 2000; Jastrebski and Robinson 2004; Millar 
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et al. 2006). Unexplained variation in these methods is often 
attributed to the vagaries of historical evolutionary events 
or unmeasured ecological variables. The role of historical 
evolutionary process on contemporary trait variation can be 
tested by evaluating the correlation between variation at neu- 
tral genetic markers and traits, thereby supporting the notion 
that historical events that influenced neutral markers also in- 
fluenced trait evolution (Merila' 1997; Merila and Crnokrak 
2001; Hankison and Ptacek 2008). A failure to detect such an 
historical correlation might be attributed to selection without 
necessarily understanding the ecological source of selection. 

More likely though, is that both contemporary and histor- 
ical processes contribute to trait variation among species and 
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populations (Taylor and McPhail 2000; Berner et al. 2010; 
Rosenblum and Harmon 201 1). Furthermore, as Langerhans 
and DeWitt (2004) point out, responses to contemporary 
selection may be influenced by historical evolutionary re- 
sponses because these shape the nature of genetic variation 
and covariation within contemporary populations. Thus, se- 
lection imposed by similar contemporary ecological condi- 
tions may result in different evolutionary responses depend- 
ing on evolutionary history. Under this hypothesis, contem- 
porary adaptive phenotypic evolution is dynamically influ- 
enced by evolutionary history. 

Disentangling the effects of historical from contemporary 
processes and the dynamic effects of history interacting with 
contemporary processes is possible when a phylogeny accu- 
rately describes the evolutionary relationships among species; 
or within species, population genetics is used to establish 
patterns of genetic structure among populations. Significant 
historical processes are expected to cause phenotypic sim- 
ilarities among populations nested in lineages that reflect 
patterns of common ancestry (Felsenstein 1985; Harvey and 
Pagel 1991; Langerhans and DeWitt 2004; Riopel et al. 2008). 
A significant correlation between population traits and en- 
vironmental features supports the role of selection imposed 
by contemporary ecological conditions; however, in the ab- 
sence of phylogenetic context, such a pattern might be var- 
iously consistent with historical or contemporary processes. 
On the one hand, a dominant role for evolutionary history 
is supported when phenotypically divergent ecomorphs oc- 
cupy similar habitats within a species and ecomorphs are de- 
scended from different ancestral lineages (Bernatchez 1997). 
In this case, there is limited evidence for the role of adaptation 
to contemporary ecological conditions since the association 
between habitat and phenotype is unreplicated, and adap- 
tive arguments must be based solely on logic derived from 
biophysical first principles. On the other hand, the parallel 
evolution of similar traits in similar environments among 
populations representing independent evolutionary lineages 
suggests a dominant role for adaptation to contemporary 
conditions (Taylor and Bentzen 1993; Bernatchez et al. 1996; 
Pigeon et al. 1997; Gislason et al. 1999; Taylor and McPhail 
1999; Schluter et al. 2004; Ostbye et al. 2006). Interactions be- 
tween historical and contemporary processes can be detected 
by statistical models that explicitly evaluate if trait variation 
among populations is affected by ancestry (using population 
genetics), contemporary environmental conditions, and in- 
teractions between ancestry and environment (Langerhans 
and DeWitt 2004). 

Studies investigating the phylogeography (Wilson and 
Hebert 1996; Bernatchez and Wilson 1998; Poissant et al. 
2005; Aldenhoven et al. 2010), ecology (Mandrak and Cross- 
man 1992; Robinson et al. 2000; McCairs and Fox 2004; Par- 
sons and Robinson 2007), and phenotypic evolution (Pigeon 
et al. 1997; Gislason et al. 1999; Jastrebski and Robinson 2004) 



of fish inhabiting postglacial lakes in the northern hemi- 
sphere have provided many valuable insights regarding how 
historical and contemporary processes interact to influence 
patterns of genetic and phenotypic variation among popu- 
lations. However, few studies have explicitly combined pop- 
ulation genetic and phenotypic variables in a single analysis 
that permits a quantitative assessment of the relative impor- 
tance of historical and contemporary factors on phenotypic 
variation among diverging populations (Ward and Mclennan 
2009; Berner et al. 2010). Here, we evaluate how contempo- 
rary and historical evolutionary processes interact to influ- 
ence the intraspecific phenotypic diversity of pumpkinseed 
sunfish (Lepomis gibbosus) in a postglacial landscape. 

Pumpkinseed sunfish polyphenism 

The pumpkinseed sunfish of North America exemplifies how 
a variety of mechanisms contribute to phenotypic diversity. 
Pumpkinseed sunfish are a common species with a wide 
distribution in northeastern North America that is more 
northerly than all other sunfish, including a substantial re- 
gion of postglacial environment (Scott and Crossman 1998). 
Its functional morphology includes specialized traits for feed- 
ing on large armored invertebrates such as snails (Lauder 
1983; Wainwright et al. 1991), which likely provides a com- 
petitive refuge from other sunfish (Mittelbach 1984). In the 
absence of bluegill sunfish (L. macrochirrus), a zooplank- 
ton generalist, some pumpkinseed sunfish are found in open 
water (pelagic) habitats of postglacial lakes, where they feed 
on zooplankton and are phenotypically divergent from local 
pumpkinseed sunfish that inhabit the shallow inshore littoral 
habitat (Robinson et al. 1993; Robinson et al. 2000; Gillespie 
and Fox 2003; Jastrebski and Robinson 2004). Phenotypic 
diversity within numerous lakes in New York State (USA) 
and the province of Ontario (Canada) is thus represented by 
a trophically related polyphenism bounded by sympatric lit- 
toral and pelagic ecomorphs (Robinson et al. 1993; Robinson 
et al. 2000; Gillespie and Fox 2003; Jastrebski and Robinson 
2004). 

We have been exploring the mechanisms that shape phe- 
notypic diversity in pumpkinseed sunfish that inhabit post- 
glacial lakes. Despite physical and ecological differences 
among lakes in this study system (Robinson et al. 2000), con- 
sistent associations between morphology (body shape and 
gill-raker length) and habitat (Robinson et al. 2000; Gille- 
spie and Fox 2003; Jastrebski and Robinson 2004) and trade- 
offs between foraging ability and phenotype (Robinson et 
al. 1996) both suggest that within-lake patterns of diversity 
are influenced by contemporary patterns of diversifying selec- 
tion between littoral and pelagic habitats that are qualitatively 
similar among lakes. The influence of historical evolutionary 
processes on shaping the pumpkinseed polyphenism though, 
has not been formally considered. 
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Historical evolutionary processes may also shape pheno- 
typic variation here because polyphenism differs between 
related species of bluegill and pumpkinseed sunfish (Ri- 
opel et al. 2008) and by meristic variation among pump- 
kinseed sunfish over a relatively large geographic scale (west- 
ern populations have higher numbers of anal and dorsal fin 
rays than eastern populations; Scott and Crossman 1998). 
This pattern is potentially explained by the historical iso- 
lation of ancestral populations in multiple glacial refugia 
(Pielou 1991). Pumpkinseed sunfish were likely present in 
both the Mississippian and the Atlantic refugia during the 
Wisconsian glaciation event (Underhill 1986; Mandrak and 
Crossman 1992) and both central Ontario and New York 
state could have been colonized by ancestral populations 
from either, or both, of the Atlantic or Misissippian refu- 
gia (Smith 1985; Underhill 1986; Mandrak and Crossman 
1992). Thus, the allopatric divergence of ancestral pump- 
kinseed populations that subsequently colonized postglacial 
aquatic habitats throughout this region over the last 15,000 
years could account for some of the within lake diversity. 
Despite an appreciation of diverse ultimate and proximate 
mechanisms shaping phenotypic diversity in pumpkinseed 
sunfish, we know little about the potential role of any histor- 
ical processes in shaping phenotypic variation among pop- 
ulations. If ancestral pumpkinseed diverged and contempo- 
rary lake populations were colonized by different or multi- 
ple ancestral sources, then some of the phenotypic diversity 
currently expressed among- and/or within-lake populations 
may reflect these historical, in addition to contemporary, 
evolutionary processes. The purpose of this study is to test 
for such historical effects using a combination of popula- 
tion genetic and morphological methods. Specifically, we 
evaluate two hypotheses about the role of historical effects. 

(1) An extreme possibility is that the phenotypic divergence 
between littoral and pelagic ecomorphs occurred only once in 
allopatry (potentially between glacial refugia), and that de- 



scendants of these populations dispersed and subsequently 
colonized lakes resulting in contemporary populations of 
polyphenic pumpkinseed sunfish. This single historical ori- 
gin hypothesis predicts that pumpkinseeds sampled from 
littoral and pelagic habitats from the same lake should be 
distantly related, and that pumpkinseeds sampled from the 
same habitat type from different lakes should be closely re- 
lated. We evaluate these predictions using allelic variation at 
six microsatellite loci. 

(2) A less extreme possibility is that ecomorphs within lakes 
have multiple contemporary origins, but that more recent 
historical effects nonetheless influence broad-scale patterns 
of phenotypic variation across the postglacial geographical 
landscape. This multiple historical origin hypothesis predicts 
that populations of pumpkinseed with different evolutionary 
histories have different phenotypes, regardless of contempo- 
rary mechanisms that currently structure diversity within 
lakes. Additionally, historical effects and adaptation to con- 
temporary habitats do not preclude each other, and might 
interact (Langerhans and DeWitt 2004). We use variation in 
allelic frequencies at six microsatellite loci to group popula- 
tions into genetic clusters with shared evolutionary histories. 
We then test for the effects of contemporary selection (sim- 
ilarities in within-lake phenotypic divergence among clus- 
ters), history (phenotypic differences between genetic clus- 
ters), and interactions between history and selection (unique 
aspects of within-lake divergence among clusters) on pheno- 
typic diversity (Langerhans and DeWitt 2004). 

Methods 
Field collections 

Pumpkinseed sunfish were collected between May and 
September 2002 from 12 lakes (hereafter referred to as "lake 
populations"), in which pumpkinseeds inhabit both the lit- 
toral and pelagic habitats (Table 1; Fig. 1). We created a 



Table 1. Locality information and sample sizes for the pelagic and littoral subpopulations inhabiting the 12 study lakes. Numbers in parentheses 
beside sample sizes (n) are the numbers of individuals in the subsample used for phenotypic analyses (see Methods). 

Lake (Abbv.) Longitude Latitude Region Drainage n pelagic n littoral 



Ashby (AS) 


45'05'N 


77°21'W 


Ontario 


Ottawa 


48(19) 


48(21) 


Mayo (MA) 


45'02'N 


77°35'W 


Ontario 


Ottawa 


48(20) 


48(18) 


Salmon Trout (ST) 


45°11'N 


77"49'W 


Ontario 


Ottawa 


48(20) 


48(23) 


Looncall (LC) 


44°45'N 


78"11'W 


Ontario 


Lake Ontario 


48(19) 


48(21) 


Monck (MO) 


45°00'N 


78"07'W 


Ontario 


Lake Ontario 


48(20) 


34(20) 


Shadow (SH) 


44°43'N 


78°48'W 


Ontario 


Lake Ontario 


48(20) 


48(20) 


Rollins (RL) 


44°38'N 


74°23'W 


New York 


St. Lawrence 


47(20) 


47(21) 


Round (RU) 


44"05'N 


74"34'W 


New York 


St. Lawrence 


48(21) 


47(20) 


Rondaxe (RA) 


43-45'N 


74°54'W 


New York 


St. Lawrence 


46(20) 


48(22) 


Paradox (PA) 


43-30'N 


73°41 'W 


New York 


Hudson 


48(20) 


48(20) 


Harris (HA) 


43°58'N 


74°07'W 


New York 


Hudson 


45(22) 


48(20) 


Lewey (LE) 


43°38'N 


74°23'W 


New York 


Hudson 


48(20) 


48(21) 
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Figure 1. Schematic illustration of locations and contemporary drainage connections among the 12 lakes used in this study. In Ontario, Ashby (AS), 
Mayo (MA), and Salmon Trout (ST) drain into the Ottawa River through the Madawaska River (Ottawa R. drainage). Looncall (LC), Shadow (SH), and 
Monck (MO) drain into the north side of Lake Ontario through the Trent River (Trent R. drainage). In the Adirondack region of northern New York 
State, three lakes drain north into the south-shore St. Lawrence River drainage. Lake Rondaxe (RA) drains through the Black river into the southeast 
side of Lake Ontario; Round Lake (RU) and Rollins Pond (RL) drain into the St. Lawrence River through the Raquette River and Lake Champlain, 
respectively. Lakes that drain south through the Adirondack region empty through the Hudson River drainage to the Atlantic Ocean and include Harris 
(HA), Paradox (PA), and Lewey (LE). 



4- 



balanced sampling design by choosing three lakes within each 
of four contemporary drainage systems, two drainages in each 
of two regions (Ontario and New York, approximately mid- 
way between glacial refugia in the Mississippi drainage, and 
along the Atlantic seaboard) . Approximately 50 individuals 
from each littoral and pelagic habitat in every lake were col- 
lected using sampling techniques consistent with previous 
studies (Robinson et al. 2000; Gillespie and Fox 2003; Jas- 
trebski and Robinson 2004). Fish were collected by angling, 
or captured in fish traps (35-cm diameter x 90-cm long; 
8-cm diameter aperture). Pelagic sunfish were designated as 
those collected from habitat adjacent to deepwater, such as 
islands and shoals in the main lake basin, and from points 
that extend into the main lake basin (>2 m depth). Littoral 
sunfish were designated as those sampled from shallow wa- 
ter (<2 m depth) close to shore, and at the back of bays, 
both usually characterized by extensive macrophyte growth. 
A study of individual movement patterns revealed that both 



ecomorphs exhibit strong habitat fidelity supporting the eco- 
logical distinctness of ecomorph populations (McCairns and 
Fox 2004). The sample of fish collected from a particular 
habitat in a lake represents the ecomorph population. Thus, 
two ecomorph populations were sampled in each lake, for 
a total of 24 ecomorph populations. A section of a pectoral 
fin from every individual was removed and stored in 95% 
ethanol for genetic analysis. Approximately the first 20 fish 
collected from each habitat in each lake were euthanized us- 
ing 250 ppm clove oil, preserved in 10% buffered formalin for 
3 months then rinsed with water, and stored in 70% ethanol. 
The whole specimens were subsequently used for morpho- 
logical analysis (see next). 

Microsatellite genotyping 

DNA was isolated from pectoral fin tissue using 
a phenol/chloroform/isoamyl alcohol extraction method 
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Table 2. . Six microsatellite loci used to study population genetic structure of pumpkinseed sunfish in central Ontario and the Adirondack region 
of New York State. Included are repeat motif, annealing temperature (T A ), source species in genus Lepomis from which these loci were isolated, and 
original reference. 



Microsatellite 


Repeat motif 


T A CC) 


Source species 


Reference 


Lma 29 


(GT)„ 


S8 


L. macrochirus 


Colbourne et al. (1996) 


RB7 


(GATA)„ 


58 


L. auritus 


DeWoody etal. (1998) 


RB20 


(GATA)„ 


58 


L. auritus 


DeWoody etal. (1998) 


Lmar9 


(AGAT)„ 


60 


L. marginatus 


Schable et al. (2002) 


Lmar14 


(AGAT)„ 


62 


L. marginatus 


Schable etal. (2002) 


Lmar18 


(AGAT)„ 


58 


L. marginatus 


Schable et al. (2002) 



(Bardacki and Skibinski 1994). Six loci with reproducible 
and clear amplification profiles (Table 2) were selected and 
amplified using the polymerase chain reaction (PCR) in 10.7 
/xl reactions (5 fi\ DNA [6 ng//xl)] with 5.7/xl PCR cocktail). 
The PCR cocktail consisted of 10 x PCR buffer (20 mM Tris- 
HC1, 50 mM KCL), 1.5 mM MgCl 2 , 0. 136 fiM dNTP, 0.36 iiM 
dCTP, 0.45 /xM of forward and reverse primers, 0.09 mg/mL 
BSA, 0.25 units of Taq DNA polymerase. Three loci (Lma 29, 
RB20, and RB7) were amplified under the following temper- 
ature regime: 95°C for 5 min; 35 cycles of 95°C for 30 sec, 
annealing temperature (T A ) for 1 min, 72°C for 1 min; then 
another 95°C for 30 sec, T A for 1 min, and finally 72°C for 
10 min. The remaining loci (Lmar 9, Lmarl4, and Lmarl8) 
were amplified using a shorter program: 95°C for 5 min; 35 
cycles of 95°C for 30 sec, T A for 45 sec, 72°C for 1 min. The 
final cycle was followed by 72° C for 5 min. The PCR products 
were separated on 6% acrylamide gels in TBE buffer and then 
visualized with a Hitachi FMBIO fluorescent imaging system 
(Hitachi, San Francisco, CA, USA) Allele sizes were deter- 
mined, using FMBIO software, by comparing the fragment 
size from the samples to several known lane standards (350 
or 500 TAMRA depending on allele size) included with every 
gel. 

Microsatellite polymorphism and 
within-lake genetic differentiation 

MSA software (Dieringer and Schlotterer 2002) was used to 
calculate allele frequencies, number of alleles, and observed 
and expected heterozygosity for each locus in every ecomorph 
population. Deviation from Hardy- Weinberg proportions 
within each ecomorph population was tested by the exact 
test for each locus and population using GENEPOP version 
3.3 (Raymond and Rousset 1995). The same program was 
used to test for linkage disequilibrium between loci within 
ecomorph populations, and to implement Fisher's exact test 
of allelic frequency homogeneity between every combination 
of ecomorph populations. To specifically determine if sym- 
patric ecomorph populations are genetically heterogeneous, 
F st values were calculated between ecomorphs in each lake 
using FSTAT (Goudet 1995). Population differentiation was 



considered significant if the confidence interval for the multi- 
locus estimate of F st , resulting from 1000 data permutations, 
excluded zero. To compare allelic richness between fish from 
different regions and drainages, f-tests (with Bonferroni ad- 
justed P-values) were used to compare the number of alleles 
in each population, averaged over all six loci. 

Population genetic structure 

There was very limited evidence of genetic divergence be- 
tween ecomorphs within lakes (see Results); therefore, eco- 
morph populations from particular lakes were pooled for all 
following analyses of geographic patterns of genetic struc- 
ture. Patterns of genetic structuring were assessed using both 
f-statistics, which do not incorporate differences in allele 
size, and ^-statistics, which incorporate the number of re- 
peated units in a microsatellite as information (assuming 
the stepwise model of microsatellite mutation). We used Ar- 
lequin (Schneider et al. 2000) to calculate pairwise values of 
F st and R st between all populations, significance was assessed 
with 1000 permutations. The allele size permutation test was 
implemented in the program SPAGeDi 1.1 (Hardy and Veke- 
mans 2002). This procedure performs 2000 permutations of 
a given dataset with randomized allele sizes, and then com- 
pares observed R st values to the distribution of permutated 
values, ^-statistics are appropriate when observed Global R st 
is improbably large (P < 0.05) compared to the permutated 
R st values; F -statistics are appropriate when the observed R st 
value is similar to the permutated R st values. 

We tested the effect of contemporary drainage distances 
on geographical patterns of among lake population genetic 
structure. First, to partition allelic variation among different 
geographic effects, we performed an analysis of molecular 
variance (AM OVA) implemented in ARLEQUIN using both 
F st and R st values. For this AMOVA, lake populations (eco- 
morphs pooled) were nested within watersheds. Second, we 
tested for isolation by distance (IBD) by evaluating the rela- 
tionship between pairwise genetic distance and geographical 
distance between populations. Google Earth was used to cal- 
culate fluvial distance among all lakes in the St. Lawrence, 
Ottawa, and Lake Ontario drainages. The lake populations 
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from the Hudson River drainage were not included in the 
IBD analyses, because contemporary gene flow between the 
Hudson populations and the rest of the study lakes along 
contemporary connections would require the movement of 
pumpkinseed through the marine environment (pumpkin- 
seed are not known to inhabit saltwater). Pairwise F st and R st 
values (see above) were transformed to F st /1-F st and R st /1-R s t 
(Rousset 1997) and the program IBDWS (Jensen et al. 2005) 
was used to perform separate Mantel tests for each distance 
measure (10,000 permutations) to determine if there was a 
significant IBD relationship. 

Population clusters were assessed using two separate ap- 
proaches. (1) We used Cavalli-Sforza and Edward's (1967) 
chord distance (D ce ) to construct a consensus neighbor- 
joining phenogram among lake populations using software 
contained in the PHYLIP 3.5 package (Felsenstein 1993). 
Statistical support for nodes was obtained from 1000 boot- 
strapped replicates of the original dataset. When reconstruct- 
ing historical population relationships between recently di- 
verged populations, studies have shown that D ce , which does 
not assume equal population sizes or mutation rates between 
loci, is more likely to achieve the correct tree topology com- 
pared to other methods (Takezaki and Nei 1996; Goldstein 
and Pollock 1997). Thus, we use D ce because we are chiefly 
interested identifying the most accurate branching topology 
for our phenogram, so that groups of populations sharing 
a recent common ancestor can be identified for morpho- 
logical analyses (see next). (2) We also used the program 
STRUCTURE 2.1 (Pritchard et al. 2000) to infer the num- 
ber of genetic clusters independently of geographic sampling 
using a posteriori Bayesian genotype clustering. We used the 
admixture model with a burn-in period of 10,000 and 10,000 
iterations and tested models with values of K (number of 
clusters) ranging from 1 to 20. The correction of Evanno et 
al. (2005) was implemented to determine the most probable 
value of K from the posterior probabilities for each K. We 
also assessed contemporary geographical population struc- 
ture by visualizing the results of STRUCTURE models with 
constrained levels of K = 2 and K = 3 (Gamier et al. 2004; 
Aldenhoven et al. 2010). 

Influence of evolutionary history on 
phenotypic variation and divergence 

The results of the above analyses of population genetic struc- 
ture informed further tests of the effect of evolutionary his- 
tory on contemporary within-lake divergence of pumpkin- 
seed ecomorphs. Specifically, we implemented the statistical 
tests suggested by Langerhans and DeWitt (2004), which in- 
volve quantifying shared and unique responses of different 
evolutionary lineages to a similar environmental gradient 
(Langerhans et al. 2003, 2006; Foster et al. 2008; Langerhans 
and Makowicz 2009; Ward and Mclennan 2009; Ruehl et al. 



20 1 1 ). In our case, the results of our population genetic analy- 
ses persistently suggest that Ontario populations of pumpkin- 
seed are distinct from New York lake populations, which are 
further subdivided by drainage — Hudson versus St. Lawrence 
(see Results) . These particular genetic clusters were useful be- 
cause they each contained replicate lakes (at least three lakes 
within each clusters) that allowed us to statistically evaluate 
the effect of history on phenotypic divergence by considering 
divergence in multiple lakes. Earlier work in this study system 
has documented parallel (shared) divergence of littoral and 
pelagic ecomorphs among lake populations of pumpkinseed 
(Robinson et al. 2000; Jastrebski and Robinson 2004). The 
objective of the current analysis is to determine if the diver- 
sity in evolutionary histories (among Ontario, St. Lawrence, 
and Hudson genetic clusters), documented in this study, is 
associated with contemporary phenotypic variation; and, to 
test for unique aspects of divergence among lake populations 
nested within the different genetic clusters. 

Variation in body shape was assessed using geometric mor- 
phometries. Photographs of each fish were taken with a Nikon 
950 digital camera (Nikon Corp., Tokyo, Japan). TPSDIG 
(Rohlf 2003) was used to locate the X and Y coordinates of 
16 homologous landmarks on the left side of each fish (Fig. 
2). TPSRELW (Rohlf 2003) was used to rotate, translate, and 
scale landmark coordinates through generalized least squares 
superimposition. From these superimposed landmark coor- 
dinates, TPSRelw also calculates affine (two uniform compo- 
nents) and nonaffine (26 partial warps) components of shape 
variation (Parsons et al. 2003). 

Nested MANCOVA models, with all partial warps and uni- 
form components as dependent variables, were used to test 
for shared and unique features of divergence (Langerhans et 
al. 2003, 2006; Langerhans and DeWitt 2004; Hendry et al. 
2006; Ward and Mclennan 2009). Throughout, we tested for 
allometry using centroid size as a covariate. First, we imple- 
mented a model that considered all 12 lakes from the three 
genetic clusters (Ontario, St. Lawrence, and Hudson). We 
tested for shared aspects of divergence (between littoral and 
pelagic ecomorphs) among lake populations using ecomorph 
(pelagic or littoral) as a fixed factor. We also tested for the 
influence of evolutionary history on phenotypic variation us- 
ing genetic cluster (Ontario, St. Lawrence, or Hudson) as a 
fixed factor. The influence of evolutionary history on diver- 
gence (unique aspects of divergence) was tested by including 
an interaction term between ecomorph and cluster. We tested 
for the influence of lake on phenotypic variation using lake 
nested within genetic cluster as a factor. The relative influ- 
ence of each factor in the model was estimated using Wilks' 
partial variance statistic {r} 2 ). Canonical variates associated 
with each effect of interest were retained and used to produce 
thin-plate spline diagrams using TPSREGR (Rohlf 2003) to 
visualize morphological variation. Note that shape deforma- 
tions along these canonical axes do not represent the original 
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Figure 2. Illustration of the 1 6 landmarks 
placed on the left side of each pumpkinseed 
sunfish, L. gibbosus, collected from 12 lakes in 
central Ontario and the Adirondack region of 
New York state, for use in the morphometric 
analysis of external body shape. 




shape space and so are cautiously interpreted (Foster et al. 
2008). The above model detected a significant effect of evo- 
lutionary history on phenotypic divergence (genetic cluster 
x ecomorph interaction — see Results) and so we performed 
additional MANCOVA models intended to investigate pat- 
terns of divergence within genetic clusters (considering lakes 
within each of the three clusters separately) . Specifically, using 
the uniform components and partial warps of morphological 
variation as dependent variables, we tested for the effects of 
allometry (using centroid size as a covariate), shared aspects 
of divergence (parallel diversity) among lakes within a cluster 
(ecomorph as a factor), the influence of evolutionary history 
on phenotypic variation (lake as a factor), and the influence 
of evolutionary history on divergence (interaction between 
ecomorph and lake as a factor). 

Results 

Microsatellite polymorphism and 
within-lake genetic differentiation 

The six microsatellite loci all showed moderate to high levels 
of polymorphism (Table SI). In many populations, H e was 
above 0.90 for certain loci. Observed heterozygosity was con- 
sistently lower than expected in the Ontario ecomorph pop- 
ulations for Lmarl4 (10/12 populations). All other loci ex- 
hibited approximately the number of deviations from HWE 
than would be predicted by chance (at a = 0.05) (5/120 pop- 
ulations). There was no significant linkage disequilibrium 
detected between pairs of loci in any ecomorph population. 
T-tests between the mean numbers of alleles per population 
averaged over all six loci indicated that the six New York pop- 
ulations had more alleles than the six Ontario populations (P 
< 0.005). No significant differences in the number of alleles 
were found between drainages after adjustment of P-values 



following the Bonferoni correction. A number of private al- 
leles shared between local ecomorph pairs (the littoral and 
pelagic ecomorphs found in the same lake) were also discov- 
ered in five Adirondack lakes: Lma29 allele 105 is only found 
in the Harris lake population; RB20 allele 303 is only found 
in the Lake Rondaxe population; Lmarl8 alleles 283 and 287 
are only found in Paradox lake; RB7 allele 314 is only found 
in Round lake; Lmar 1 4 allele 4 1 1 is only found in Lewey lake. 

The null hypothesis of homogeneity of allele frequencies 
was rejected, using Fisher's exact test, in all pair- wise compar- 
isons of ecomorph populations (P < 0.001), excluding com- 
parisons between sympatric ecomorphs within lakes. There 
was only very limited evidence of genetic differentiation be- 
tween pelagic and littoral ecomorphs within lake populations. 
Between local ecomorphs, the null hypothesis of homogene- 
ity of allele frequencies was rejected (P < 0.05) in two lakes: 
Round Lake in the St. Lawrence drainage and Lewey Lake in 
the Hudson drainage. However, multilocus pairwise F st val- 
ues between local ecomorphs were very low (<0.008), and 
only in one lake (Round Lake) was it significantly different 
from zero (F st = 0.005). Thus, ecomorph populations from 
the same lake were pooled in order to investigate broader 
geographic patterns of genetic structure. 

Population genetic structure 

Allele size permutation tests revealed that geographical vari- 
ation in allele sizes contributed useful phylogeographic in- 
formation (P < 0.05 in all tests); therefore, we present the 
results of both F st - and R st -based distances. All multilocus, 
pairwise estimates of F st and R st were greater than zero (P 
< 0.01) (Table 3). The hierarchical (AMOVA) analyses of 
genetic variance revealed that allelic variation is attributable 
to geographical entities (Table 4). We detected significant 
patterns of population structure reflecting lake and drainage 
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Table 3. . Estimates of multilocus pairwise R st (above diagonal) and F 9 (below diagonal) values based on variation at six polymorphic microsatellite 
loci among the 12 lake populations of pumpkinseed sunfish (lake abbreviations found in Table 1). All values were significantly greater than zero (P < 
0.01). 





AS 


MA 


ST 


MO 
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SH 


RL 


RU 


RA 


PA 


LE 


HA 


AS 
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0.120 


0.159 


0.094 


0.098 


0.050 


0.371 


0.792 


0.405 


0.112 


0.177 


0.254 


MA 


0.095 




0.293 


0.123 


0.220 


0.144 


0.333 


0.778 


0.369 


0.121 


0.281 


0.223 


57 


0 1 04 


0 137 




0 147 


0 111 


0 114 


0 290 


0 782 


0 311 


0 125 


0 391 


0 1 58 


MO 


0.160 


0.180 


0.095 




0.077 


0.044 


0.294 


0.781 


0.325 


0.050 


0.288 


0.167 


LC 


0.099 


0.109 


0.097 


0.144 




0.029 


0.279 


0.765 


0.319 


0.038 


0.227 


0.156 


SH 


0.040 


0.060 


0.091 


0.135 


0.061 




0.296 


0.766 


0.328 


0.040 


0.207 


0.170 


RL 


0.138 


0.109 


0.157 


0.183 


0.102 


0.093 




0.228 


0.016 


0.242 


0.457 


0.042 


RU 


0.153 


0.138 


0.172 


0.193 


0.113 


0.110 


0.025 




0.274 


0.684 


0.797 


0.445 


RA 


0.177 


0.176 


0.198 


0.217 


0.128 


0.150 


0.048 


0.033 




0.250 


0.490 


0.069 


PA 


0.133 


0.129 


0.134 


0.170 


0.068 


0.085 


0.061 


0.066 


0.080 




0.192 


0.133 


LE 


0.131 


0.129 


0.149 


0.180 


0.114 


0.094 


0.072 


0.078 


0.111 


0.083 




0.382 


HA 


0.068 


0.075 


0.077 


0.114 


0.056 


0.039 


0.048 


0.042 


0.074 


0.050 


0.060 





Table 4. . Hierarchical AMOVA based on six microsatellite loci in pop- 
ulations of pumpkinseed sunfish from lakes in central Ontario and the 
Adirondack region of New York State. AMOVA was used to test the 
hypothesis that contemporary drainage structure explains patterns of 
genetic variation by pooling ecomorphs within lakes and nesting lakes 
within drainages. For each parameter, upper values are based on allele 
frequency only (F st ) and lower values considered variation in allele size 
and allele frequency (R st ). Listed are variance (V), proportion of variance 
explained by a particular grouping (%), fixation index (haplotypic corre- 
lation at corresponding grouping), and the probability of having a more 
extreme variance component than that observed (P). 



Effect 


V 


% 


Fixation 
indices 


P-value 


Among drainages 


0.093 


3.93 


0.039 


<0.001 




40.7 


27.2 


0.037 


<0.001 


Among lakes within 


0.18 


7.61 


0.079 


<0.001 


drainages 


20.5 


13.7 


0.081 


<0.001 


Among individuals 


2.1 


88.4 


0.115 


<0.001 


within lakes 


88.5 


59.1 


0.115 


<0.001 



membership. We note that in all AMOVA tests of genetic 
structure between geographical entities (as opposed to within 
populations), calculations of R st explained more variation 
than calculations of F st . When only F st calculations were con- 
sidered, allelic differences among drainages explained less 
variation than allelic variation among lakes within drainages 
(3.9% and 7.6%, respectively); however, when only R st val- 
ues were considered, variation among drainages explained 
a greater proportion of molecular variation than variation 
among lakes within drainages (27.2% and 13.7%, respec- 
tively). Despite the AMOVA results (which detected signif- 
icant effects of geographical entities), mantel tests compar- 
ing genetic (R st /1-R st and F st /1-F st ) and geographic distance 
matrices found no significant pattern of IBD among these 
lake populations (St. Lawrence, Ottawa, and Lake Ontario 



drainages: all P > 0.1). Taken together, pumpkinseeds sam- 
pled from particular lakes and drainages seem to be closely re- 
lated to each other, however, these affiliations do not seem to 
reflect distances among contemporary waterways (see next). 

The neighbor-joining phenogram based on D ce distances 
had a robust topology with several nodes being supported in 
over 70% of the bootstrap replicates (Fig. 3 ) . All Ontario pop- 
ulations form a reasonably well-supported cluster separate 
from the New York populations (78.9% of bootstrap repli- 
cates). In New York, the populations corresponding to the 
St. Lawrence drainage form a well-supported (100% boot- 
strap support) grouping, separate from the Hudson River 
drainage populations — which were more closely related to 
Ontario populations. In contrast, populations within con- 
temporary drainage systems in Ontario do not form mono- 
phyletic groupings; for example, populations from Monck 
and Salmon Trout lakes, which are from different Ontario 
drainages, form a monophyletic cluster with 97.8% boot- 
strap support. 

Individual-based clustering using STRUCTURE generally 
identified population clusters associated with individuals 
sampled from the 12 different lakes. However, the most prob- 
able number of clusters, based on the approach of Evanno et 
al. (2005), in our dataset was actually K = 14 (Table 5), 
indicating more population genetic clusters than lakes. Un- 
der the K = 14 STRUCTURE analyses, two clusters (5 and 
12) contain relatively high numbers of individuals from dif- 
ferent lakes. Cluster 5 includes individual fish from Rollins, 
Round, and Rondaxe lakes from the Hudson River drainage, 
cluster 12 includes individual fish from Ashby and Shadow 
lakes from Ontario (Fig. 4; Table 6). We interpret these extra 
clusters to be an artifact of how recently these lakes have 
diverged from each other, and we are not aware of any 
recent introductions or transfers of pumpkinseeds among 
lakes that could produce this pattern. Furthermore, these 
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New York i Ontario 



Figure 3. Phenogram depicting relationships 
among the 12 lake populations of 
pumpkinseed sunfish sampled from Central 
Ontario and the Adirondack region of New 
York State. Superscripts correspond to 
drainage systems: Ottawa R. drainage of 
Ontario (■); Trent R. drainage of Ontario (a); 
St. Lawrence R. drainage of New York (•); 
Hudson R. drainage of New York (▼). 
Relationships were inferred from a matrix of 
D cli genetic distances and the Neighbor-joining 
algorithm. For 1000 bootstrap replicates, node 
values of 50% and higher are shown. This 
analysis grouped population from Ontario 
together with 78.9% bootstrap support, 
populations from the St. Lawrence drainage of 
New York with 100% bootstrap support. 
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Figure 4. Bar plot of membership proportions 
for all individual from the 12 lakes (represented 
by abbreviations, see methods) in this study 
using STRUCTURE for K = 2, K = 3, and K = 1 4 
(most likely). Symbols correspond to drainage 
systems: Ottawa R. drainage of Ontario (■); 
Trent R. drainage of Ontario (a); St. Lawrence R. 
drainage of New York (•); Hudson R. drainage of 
NewYork (t). 
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lakes were not those with marginal support for allelic het- 
erogeneity between ecomorphs, suggesting that these genetic 
clusters do not reflect differences in ecological specializa- 
tion. Individuals from all other lakes in the K = 14 model 



were generally assigned to lake-specific clusters with rela- 
tively high frequency (Table 6). We also used STRUCTURE 
to assess whether the results of the Bayesian individual-based 
clustering approach generally corroborated the more basal 
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Table 5. . Likelihood of various numbers of genetic clusters (K) based on 
the STRUCTURE 2.2 analysis of all 1 129 individual pumpkinseed sunfish 
from the 12 lake populations considered in this study. Following Bayes' 
rule the probability (PrK) for K = 14 was estimated to be 1 . 



K 


In Pr (X IK) 


PrK 


1 


-31,640.7 


0 


2 


-29,765 


0 


3 


-28,887.1 


0 


4 


-28,299.5 


0 


5 


-27,777.9 


0 


6 


-28,633.8 


0 


7 


-27,254 


0 
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-27,045.4 
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-26,963.2 
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1 0 


-26 782 3 
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1 1 


-26,719.2 
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12 


-26,642 
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13 


-26,790.5 


0 


14 


-26,528.8 


1 


15 


-26,764.5 


0 


16 


-29,727.2 


0 


17 


-26,866.9 


0 


18 


-300,066.1 


0 


19 


-27,138.1 


0 


20 


-28.410.2 


0 



geographical groupings suggested by the D ce phenogram. 
Specifically, when we constrained STRUCTURE to estimate 
two clusters (K = 2), individuals from Ontario were generally 
distinguished from individuals from New York (Fig. 4), con- 
sistent with the D ce phenogram. Also, when we constrained 
STRUCTURE to estimate three clusters (X = 3), individu- 
als from New York were further subdivided into drainage- 
specific clusters (Fig. 4), consistent with the D ce phenogram. 

Distributions of allelic frequencies suggest that New York 
populations are distinct from Ontario populations, and that 



lake populations within the St. Lawrence drainage within 
New York are closely related to each other (histograms of 
allelic frequencies in Fig. SI). For all loci, the allele size range 
present in Ontario was present in New York. However, New 
York populations often contained an additional size range not 
found in Ontario. The allelic size range appeared bimodal for 
several loci: Lma29, Lmarl8, RB7, and Lmarl4. Geographic 
variation in allele size range was most dramatic for RB7. In 
this case, an allele size range of approximately 150-200 base 
pairs characterized Ontario populations, while an additional 
mode of alleles (250-300 base pairs) was common in New 
York populations. 

Shared and unique aspects of within-lake 
phenotypic divergence 

The full MANCOVA analysis of all 12 lakes from the three 
genetic clusters detected a significant effect of ecomorph on 
morphological variation (Table 7), indicating that there is a 
shared aspect of divergence between littoral and pelagic habi- 
tats among all lakes. Pelagic ecomorphs had generally more 
fusiform bodies, narrower pectoral fin insertions, and per- 
haps larger eyes and mouth regions than littoral ecomorphs 
(Fig. 5). This result generally confirms the findings of ear- 
lier work in this system. However, our findings also demon- 
strate substantial additional variation in external body form 
in pumpkinseed sunfish that can be related to evolutionary 
history. 

In the full analysis of all 12 populations (Table 7) and apart 
from allometry (partial variance centroid size = 72.2%), the 
effect of genetic cluster accounted for the largest propor- 
tion of phenotypic variation (partial variance genetic cluster 
= 60.4%), a much higher proportion of phenotypic varia- 
tion than was accounted for by the shared response to con- 
temporary divergent selection between lake habitats (partial 



Table 6. . Results from STRUCTURE based on the most likely number of cluster (K = 14). Proportion of genotypes from each sampling site within 
each of the 14 inferred clusters are shown. Bold values indicate proportions greater than 0.1 . 
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Lake 1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


1 1 


12 
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AS 


0.017 


0.009 
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0.014 


0.007 


0.019 


0.031 


0.02 
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0.089 


0.318 


0.025 
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0.019 


0.014 
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0.011 


0.019 


0.028 
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0.03 
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0.02 
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0.009 


ST 
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0.018 


0.014 


0.006 
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0.016 


0.01 


0.035 


0.023 


0.052 
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0.007 


MO 
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0.006 


0.011 
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0.006 
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0.012 
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0.017 
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0.034 


0.006 
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0.012 


0.034 


0.021 


0.012 
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0.049 
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0.053 
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SH 
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0.05 
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0.086 
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0.017 
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0.028 
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Table 7. . MANCOVA analyses examining how shape variation is related to various factors. As a result of finding significant interactions between 
ecomorph and genetic cluster, separate MANCOVA analyses were performed for each genetic cluster (St. Lawrence, Hudson, and Ontario). Note that 
partial variance explained does not sum to unity (Langerhans and DeWitt 2004). 

df (numerator, Wilk's Partial variance 

F denominator) P partial X explained (%) 



All Populations 



Centroid size 


41.33 


28,445 


0.0001 


0.278 


72.2 


Ecomorph 


4.53 


28,445 


< 0.0001 


0.778 


22.2 


Genetic cluster 


24.28 


56,890 


< 0.0001 


0.156 


60.4 


Ecomorph x cluster 


2.58 


56,890 


< 0.0001 


0.740 


13.9 


Lake (Cluster) 


7.45 


252,3845 


< 0.0001 


0.033 


26.7 


:. Lawrence populations 












Centroid size 


8.87 


28,90 


< 0.0001 


0.266 


73.4 


Ecomorph 


3.24 


28,90 


< 0.0001 


0.498 


50.2 


Lake 


7.50 


56,180 


< 0.0001 


0.090 


70.0 


Ecomorph x lake 


1.48 


56,180 


< 0.0283 


0.469 


31.5 


udson populations 












Centroid size 


8.11 


28,89 


< 0.0001 


0.110 


89.0 


Ecomorph 


2.17 


28,89 


0.0032 


0.594 


40.6 


Lake 


6.06 


56,178 


< 0.0001 


0.118 


65.6 


Ecomorph x lake 


1.18 


56,178 


0.214 


0.533 


27.0 


ntario populations 












Centroid size 


13.89 


28,201 


< 0.0001 


0.341 


65.9 


Ecomorph 


6.68 


28,201 


< 0.0001 


0.694 


30.6 


Lake 


7.50 


140,997 


< 0.0001 


0.029 


83.0 


Ecomorph x lake 


3.79 


140,997 


< 0.0001 


0.122 


65.1 



variance ecomorph = 22.2%), or by similarities among 
pumpkinseeds from the same lake (partial variance lake = 
26.7%). In Figure 6 we visualize the canonical variates corre- 
sponding to the effect of genetic (historical) cluster. The first 
canonical variate ( CV 1 ) indicates that pumpkinseed from the 
Hudson genetic cluster have deeper bodies, wider pectoral in- 
sertion widths, and possibly smaller mouths compared to fish 
from the St. Lawrence and Ontario genetic clusters. Variation 
along CV2 indicates that pumpkinseed from the St. Lawrence 
cluster have narrower bodies, more ventrally placed pectoral 
fins, and larger mouths than pumpkinseeds from the Hud- 
son and Ontario clusters. This is consistent with the signif- 
icant interaction between ecomorph and cluster from the 
full MANCOVA model of all 12 populations that indicates 
unique aspects of morphological divergence among the three 
genetic clusters (Table 7). Therefore, in this system, patterns 
of morphological variation are influenced by both contem- 
porary processes acting between divergent environments and 
evolutionary history. 

Given this result, we were motivated to visualize axes of 
phenotypic divergence associated with habitat for popula- 
tions in each genetic cluster, and so performed additional 
MANCOVA analyses only considering populations within 
each genetic cluster. In the St. Lawrence populations, vari- 
ation in body depth was not associated with different lake 
habitats; but instead, divergence between pelagic and littoral 
habitats involved variation in the width and position of the 



pectoral fin insertions (Fig. 5). In contrast, among the Hud- 
son and Ontario populations, the width and position of the 
pectoral fin insertion were similar for both ecomorphs (Fig. 
5), but the body shape of the Littoral ecomorph was rounder 
and more gibbose — this result is more pronounced in the 
Hudson cluster compared to the Ontario cluster. 

Discussion 

A central challenge in evolutionary biology is to understand 
how the interplay between historic and contemporary evolu- 
tionary process affects patterns of diversity within a species. 
Our population genetic analyses provide considerable in- 
sight into this interplay, because they revealed that while 
there is strong indirect evidence that diversifying selection 
acts on sympatric ecomorphs within lakes, contemporary 
lake populations of pumpkinseed also show obvious regional 
patterns of genetic structure. This raises the possibility that 
postglacial aspects of ancestry shapes phenotypic variation 
in two interesting ways. First, our morphometric analyses re- 
vealed strong associations between variation in body shape 
and genetic cluster membership. Second, unique aspects of 
within-lake phenotypic divergence occurred among genetic 
clusters, suggesting that diversification between lake habi- 
tats involves an interaction between postglacial history and 
contemporary natural selection within lakes. We discuss the 
evidence for various scales of historical divergence before 
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Figure 5. Shared features of body shape 
divergence between pelagic and littoral 
ecomorphs of pumpkinseed. Circles (filled = 
littoral, open = pelagic) indicate a particular 
population's mean canonical value extracted 
from the "ecomorph" canonical axis of the 
MANCOVA analyses (see Methods section). 
Thin-plate spline transformations represent the 
most extreme deviation from the consensus 
configuration (magnified x2), and were 
generated by regressing canonical scores 
against shape data (X-Y coordinates) using 
TPSregr. As a result of finding significant 
interactions between ecomorph and genetic 
cluster in our analysis that included all 
populations (see Results), we also visualized 
divergence within each genetic cluster (St. 




Lawrence, Hudson, and Ontario). -0.225 -0.150 -0.075 0 000 0 075 0150 0.225 



focusing on how these have affected contemporary patterns 
of phenotypic diversity in pumpkinseed sunfish. 

Within- and among-lake patterns of 
population genetic structure 

A major goal was to evaluate if contemporary populations of 
polyphenic pumpkinseed could have resulted from a single 
allopatric divergence between glacial refugia (e.g., Bernatchez 
1997). We found no evidence of multiple colonizations fol- 
lowing a single historical divergence, and instead our results 
are more consistent with within-lake divergence following 
colonization. The null hypothesis of allelic homogeneity be- 
tween sympatric ecomorph populations was not rejected in 



nearly all lakes. Round and Lewey lakes showed some evi- 
dence of allelic heterogeneity, but multilocus F st values be- 
tween ecomorphs were very low ( <0.008), and in only Round 
Lake was this different from zero (F st = 0.005). Expanded 
within-lake analyses of individual movement between habi- 
tats (McCairns and Fox 2004) may help us better understand 
the minimal neutral genetic divergence of sympatric eco- 
morphs (e.g., Fraser and Bernatchez 2005). Our individual- 
based clustering analyses using STRUCTURE also found little 
evidence of close relatedness among ecomorphs sharing the 
same lake habitat. In addition, we find private alleles shared 
between local ecomorph pairs in five populations (Harris, 
Rondaxe, Paradox, Round, and Lewey lakes) . These data cause 
us to reject the ideas that ecologically and phenotypically 
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Figure 6. Features of body shape variation 
associated with genetic cluster membership 
(Ontario, St. Lawrence, and Hudson — see 
Results). Shapes represent different genetic 
clusters (filled = littoral ecomorphs, open = 
pelagic ecomorphs), triangles are Ontario 
populations, circles are Hudson populations, 
and squares are St. Lawrence populations. 
The x -axis represents mean body shape 
variation depicted by canonical variate 1 
(CV1), they -axis represents mean body 
shape variation depicted by canonical variate 
2 (CV2), which were both extracted from 
the MANCOVA analysis that considered all 
populations (see Methods). Thin-plate spline 
transformations represent the most extreme 
deviation from the consensus configuration 
(magnified x2) and were generated by 
regressing canonical scores against shape 
data (X-Y coordinates) using TPSregr. 
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divergent pumpkinseeds within lakes are distantly related 
and also reject that similar ecomorphs from different popu- 
lations are closely related as expected under a single historical 
divergence hypothesis. 

This conclusion is provisional, however, because we can- 
not reject all possible hypotheses regarding single evolution- 
ary origins of ecomorphs within lakes without using a more 
powerful population genomic approach. For example, eco- 
morphs within lakes could arise from a single allopatric di- 
vergence followed by extensive within-lake introgression ex- 
cept at ecologically important genomic islands (Ellison et 
al. 201 1), which, under persistent divergent selection, are not 
homogenized between ecomorphs. This seems less likely here 
though because we see little, if any, divergence in neutral alle- 
les between ecomorphs in any lake populations despite a ro- 
bust dataset involving strong replication within and among 
lakes. Moreover, diversifying selection arising between lake 
habitats is a key component of this hypothesis in order to 
maintain divergent phenotypes, and so even if this process 
had occurred, a role for diversifying selection on sunfish be- 
tween lake habitats would still be supported. 

While we found little evidence that contemporary pump- 
kinseed sunfish lake ecomorphs derive from a single historical 
divergence, we nonetheless found some evidence that an- 
cestral pumpkinseed sunfish likely inhabited multiple glacial 
refugia. Complicated patterns of genetic structure and phylo- 
geography are common for North American fish inhabiting 
formerly glaciated regions. Many fish colonized large geo- 
graphical areas from various glacial refugia ( Pielou 1991; Wil- 
son and Hebert 1996; Bernatchez and Wilson 1998; Poissant 
et al. 2005; Aldenhoven et al. 2010). Specific dispersal out- 
lets produced extensive zones of secondary contact between 



refugia lineages (Bernatchez and Dodson 1991; Bernatchez 
1997; Wilson and Hebert 1996; Turgeon and Bernatchez 
2001). After secondary contact, some reproductively iso- 
lated glacial lineages persisted in sympatry (Bernatchez 1997) , 
while others resulted in extensive introgression (Turgeon and 
Bernatchez 2001; Lu et al. 2001). During the Wisconsian 
glaciation event, pumpkinseed sunfish were likely present in 
both the Mississippian and the Atlantic refugia (Underhill 
1986; Mandrak and Crossman 1992) and so may have colo- 
nized central Ontario and New York state from either or both 
(Smith 1985; Underhill 1986; Mandrak and Crossman 1992). 
Our analyses suggest that Ontario populations of pumpkin- 
seed received a greater genetic contribution from Mississip- 
pian ancestors, while New York populations received a greater 
contribution from Atlantic refugial ancestors, although in- 
trogression between lineages likely occurred in both regions. 
Pumpkinseed populations in New York contain an additional 
allele size range not present in Ontario populations that was 
also bimodal for four loci in the New York populations (Fig. 
SI). 

Our postglacial dispersal hypothesis of pumpkinseed sun- 
fish is similar to one proposed for lake trout in eastern North 
America. Wilson and Hebert (1996) found that Atlantic hap- 
lotypes dominated in New York populations while Missis- 
sippian haplotypes dominated central Ontario populations 
of lake trout. While generally classified as a warm water fish, 
pumpkinseed sunfish are tolerant of cooler waters, having the 
lowest critical thermal maxima of Lepomid taxa (Beitinger et 
al. 2000, Coker et al. 2001), and so perhaps dispersed some- 
what like the cold-tolerant lake trout across the postglacial 
landscape. A more genomically and geographically extensive 
sampling effort ideally incorporating a more nuanced under- 
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standing of periglacial waterways (i.e., Poissant et al. 2005), 
is required to further evaluate the geographical origins and 
dispersal of pumpkinseed sunfish ancestors in the region. 

A significant insight of this study is the finding that pop- 
ulation genetic structure reflected geographical entities. Our 
AMOVA analyses indicated that both lakes and drainages ex- 
plained significant allelic variation. All pairwise estimates of 
F st and R st between lakes were significantly greater than zero 
(F st range = 0.025-0.198; R st range 0.016-0.792). Also, the 
STRUCTURE analysis identified K = 14 as the most likely 
number of clusters (two more clusters than lakes; although 
clusters were generally associated with particular lakes). Yet, 
we found no evidence that the magnitude of genetic dif- 
ferentiation among lake populations was related to distance 
between lakes along contemporary watersheds (e.g., Crispo 
and Chapman 2008). This conflict largely arises from the 
tendency of Ontario populations in the Trent and Ottawa 
drainages to be closely related, despite the relative proximity 
of the Trent drainage to the St. Lawrence drainage in New 
York (Fig. 1). The young geological ages of these contem- 
porary waterways may not have provided enough time to 
develop IBD patterns (Bernatchez and Wilson 1998). 

Our population genetic analyses also persistently found 
that lake populations clustered into three larger geograph- 
ical entities: Ontario, the south shore of the St. Lawrence 
River, and the Hudson River drainages. The neighbor-joining 
phenogram based on D ce distances grouped lake populations 
to these three geographical entities and this was corroborated 
by a series of constrained STRUCTURE models ( Gamier et al. 
2004; Crispo and Chapman 2008; Aldenhoven et al. 2010). 
A model constrained to K = 2 separated sunfish between 
Ontario and New York, while a K = 3 model distinguished 
sunfish from the regions of Ontario, St. Lawrence, and the 
Hudson River drainages. The absence of detailed evidence for 
IBD along contemporary waterways above, however, suggests 
that our study lakes were colonized along different historical 
rather than contemporary waterways (Poissant et al. 2005), 
such as mobile proglacial lakes that followed the retreating 
edge of the glaciers (e.g., Mandrak and Crossman 1992). 
While this regional population genetic structure could result 
from neutral and/or unknown adaptive evolutionary pro- 
cesses, this geographic pattern is nonetheless well supported 
and provides an opportunity to evaluate its effects on contem- 
porary phenotypic diversification of pumpkinseed sunfish. 

Contemporary and historical influences on 
phenotypic variation 

Previous studies documenting parallel aspects of phenotypic 
divergence among lakes containing polyphenic populations 
of pumpkinseed (Robinson et al. 2000; Jastrebski and Robin- 
son 2004; Gillespie and Fox 2003; Riopel et al. 2008) have as- 
sumed that lake populations are evolutionarily independent 



and that historical effects have little affect on phenotypic di- 
versity. Our analyses above show that neither of these assump- 
tions is strictly justified because all populations are related, 
although in a nested hierarchy that affords less evolutionary 
independence between lakes within regional genetic clusters 
and greater independence among lakes from different genetic 
clusters. Therefore, does incorporating our understanding of 
these historical effects help us better understand the nature 
of polyphenism in pumpkinseed sunfish? We evaluated the 
influence of history on body shape variation in this system 
by: (1) comparing the phenotypes of pumpkinseeds from 
different genetic clusters (identified by our analyses of mi- 
crosatellite allele frequencies), and (2) testing for unique and 
shared aspects of within-lake divergence among these genetic 
clusters (Langerhans and DeWitt 2004; Hendry et al. 2006; 
Ward and Mclennan 2009). 

We continue to find evidence of parallel diversification in 
some aspects of body form between lake habitats after ac- 
counting for history (a significant effect of ecomorph) that 
is consistent with earlier work (Robinson et al. 2000; Jas- 
trebski and Robinson 2004). Pelagic ecomorphs generally 
have more streamlined bodies and narrower insertions of the 
pectoral fins compared to littoral ecomorphs (Fig. 5). These 
differences are consistent with established form-function re- 
lationships between fish body shape and swimming perfor- 
mance (reviewed in Jastrebski and Robinson 2004; Riopel et 
al. 2008). A more streamlined body may permit more effi- 
cient locomotion, but also reduces maneuverability (Webb 
1984; Langerhans et al. 2003). Thus, pelagic pumpkinseeds 
that forage on aggregates of zooplankton dispersed in the 
large but structurally less complicated pelagic habitat may 
be well served by a body shape emphasizing swimming ef- 
ficiency at the expense of maneuverability. Pectoral fins are 
used by sunfish for turning and braking and also contribute 
importantly to paired fin locomotion (Drucker and Lauder 
2001, 2002), especially when precise movements are required 
in structurally complex habitats (but see Highamet al. 2005). 
Thus, littoral pumpkinseeds may be well served by wider pec- 
toral fins if they reflect an increased reliance upon precise, 
pectoral fin locomotion (Ehlinger and Wilson 1988). 

Our results suggest that the littoral-pelagic polyphenism is 
replicated because lineages of pumpkinseed sunfish become 
increasingly evolutionarily independent especially between 
biogeographic regions, while polyphenism is still found in 
many lakes. Replicated phenotypic polyphenisms support the 
role of divergent selection between littoral and pelagic habi- 
tats in many of the postglacial lakes of central Ontario and 
the Adirondack region of New York State. However, there was 
little evidence of heterogeneity in allelic frequency between 
ecomorphs in most lake populations, indicating that sym- 
patric ecomorphs have not been sufficiently reproductively 
isolated for neutral alleles to become ecologically isolated at 
this scale as in other studies (Crispo and Chapman 2008). 
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While we have found evidence that ancestral populations of 
pumpkinseed sunfish likely diverged between glacial refugia 
and subsequently colonized postglacial lakes, contemporary 
processes occurring within lakes continue to shape the phe- 
notypic diversity between coexisting sunfish ecomorphs. 

Despite numerous studies that have documented associa- 
tions between neutral genetic variation and phenotypic vari- 
ation (Merila 1997; Hendry et al. 2002; Nicholls et al. 2006; 
Ward and Mclennan 2009), we were nonetheless surprised 
that the effect of postglacial colonization history (genetic 
cluster) also seems to influence body shape, and do so more 
strongly than either the effect of lake habitat (ecomorph) 
or lake population (full analysis of all 12 lake populations 
after accounting for allometry; Table 7). We visualized the 
body shape differences associated with genetic cluster mem- 
bership, and found that the first two canonical variates both 
emphasized variation in body depth. Pumpkinseeds from 
the Hudson River cluster tend to have the greatest body 
depths, while pumpkinseed from Ontario had deeper bod- 
ies than those from the St. Lawrence genetic cluster. Evolu- 
tionary history and contemporary patters of selection may 
be confounded such that the geographic entities associated 
with our genetic clusters (Ontario, Hudson river drainage, 
St. Lawrence river drainage) differ in contemporary envi- 
ronmental features. This limits our ability to conclude that 
differences in morphology associated with genetic cluster are 
exclusively caused by historical evolutionary processes. How- 
ever, we know of no obvious contemporary ecological differ- 
ence (i.e., in aquatic community structure or physical habitat 
characteristics) of lakes among these regions or drainages that 
could produce these effects (latitudinal differences between 
Ontario and New York may also be offset by higher altitudes 
of lakes in the Adirondack region of New York relative to 
Ontario). A comprehensive limnological survey of a subset 
of lakes representing each geographical region could usefully 
address this issue. Our preliminary conclusion is that it seems 
more reasonable that regional genetic clusters reflects ances- 
tral lineages that diverged under selection, genetic drift, or 
founder events either between glacial refugia or during the 
dispersal phase that preceded the colonization of contempo- 
rary lakes (Ward and Mclennan 2009). 

Evolutionary history has also influenced the contempo- 
rary process of phenotypic divergence between littoral and 
pelagic lake habitats (an interaction between genetic cluster 
and ecomorph; Table 7). To better understand patterns of 
phenotypic divergence in each regional cluster, we visualized 
the canonical variates associated with ecomorph from three 
cluster-specific MANCOVAs. Comparisons revealed some in- 
teresting insights regarding interactions between postglacial 
history and responses to contemporary selective regimes be- 
tween lake habitats. Within the St. Lawrence cluster, differ- 
ences between pelagic and littoral pumpkinseeds stress vari- 
ation in the width of the pectoral fin insertion (consistent 



with the results of the analysis that included all populations, 
St. Lawrence littoral ecomorphs had wider pectoral fin inser- 
tions). However, in contrast to our all-population analysis, 
the St. Lawrence littoral ecomorphs also seem to be more 
streamlined than their local pelagic counterparts. Within 
the Hudson and Ontario clusters, the pectoral fin insertion 
width seems to be relatively invariant between ecomorphs, 
but pelagic ecomorphs have more streamlined body shapes 
than their littoral counterparts. Thus, depending on their 
postglacial evolutionary history, pumpkinseeds can evolve di- 
verse morphological solutions to divergent selection between 
littoral and pelagic habitats. Improvements to maneuverabil- 
ity (favored in the littoral habitat) can perhaps be made by 
modifying the width of pectoral fin insertions and/or by in- 
creasing body depth. Conclusions about the exclusive effects 
of history are again provisional because some cryptic con- 
temporary ecological differences between these geographical 
entities may be confounded with historical evolutionary ef- 
fects. Assessing contemporary patterns of selection in repli- 
cated lake populations nested within the different genetic 
clusters would address this hypothesis. Selection data (the 
relationship between traits and fitness), could be combined 
into a single selection analysis (Weese et al. 2010) in order to 
test for similarities in the pattern and strength of selection 
among lakes nested within genetic clusters. 

It is important to acknowledge that since our morpho- 
logical analyses involved wild caught fish, we cannot dis- 
criminate between genetic versus environmental variation, 
although this has been the focus of substantial previous re- 
search in this system. Rearing studies demonstrate that a 
combination of minor genetic and substantial plastic de- 
velopmental responses to dietary and predator cues shapes 
pumpkinseed sunfish phenotype including between these 
ecomorphs (Robinson and Wilson 1996). Of greater interest, 
is the finding that the plastic developmental responses them- 
selves have diverged between sympatric ecomorphs (Parsons 
and Robinson 2006, 2007; Januszkiewicz and Robinson 2007; 
Robinson et al. 2008). The near absence of genetic differen- 
tiation between sympatric ecomorphs at neutral loci here 
combined with heritable differences between sympatric eco- 
morphs raises the possibility that certain functional loci may 
adaptively diverge between ecomorphs perhaps via a "porous 
genome" phenomenon (Barton and Bengtsson 1986, Wu 
2001). Under this model, early in a divergence, only func- 
tional loci under diversifying selection diverge while neutral 
loci remain homogenized for much longer between sympatric 
populations (Via and West 2008). If occurring here, then it is 
interesting that the functional loci apparently diverging be- 
tween sympatric pumpkinseed ecomorphs influence plastic 
morphological responses to littoral versus pelagic conditions. 
The role of phenotypic plasticity in adaptive divergence is in- 
completely understood (Ghalambor et al 2007; Pfennig et al 
2010; Moczek et al 2011), and it can either enhance or con- 
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strain adaptive divergence through its effect on gene flow 
among subpopulations experiencing diversifying selection 
(Price et al 2003). Evaluating the porous genome hypoth- 
esis and the role of plastic traits here will require identifying 
nonneutral genetic markers under diversifying selection be- 
tween lake habitats and then testing their relationship to the 
divergence of plastic morphological responses between eco- 
morphs. 

Our study documents how the variable evolutionary his- 
tories of different populations diverging along a similar en- 
vironmental gradient produce both shared and unique as- 
pects of divergence. In addition, our understanding of these 
new regional effects on phenotype within pumpkinseed sun- 
fish fills a gap in an array of historical effects that range 
from deep effects identified among Lepomid taxa predat- 
ing that last glacial period (Riopel et al. 2008), more recent 
historic glacial effects on meristic characters between east- 
ern and western populations of pumpkinseed sunfish over 
their range (Scott and Crossman 1998), and the contempo- 
rary processes of diversification occurring within and among 
lakes (Robinson et al. 2000; Gillespie and Fox 2004; Jastreb- 
ski and Robinson 2004). By incorporating population genetic 
with morphometric analytic methods, we have highlighted 
how contemporary patterns of phenotypic divergence are 
better understood against a known backdrop of evolutionary 
history. 
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